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Abstract 


Micro-mechanical models for a study of nonlinear visco-elastic response of composite laminae 
are developed and their performance compared. A single integral constitutive law proposed by 
Schapery and subsequently generalized to multi-axial states of stress is utilized in the study for 
the matrix material. This is used in conjunction with a computationally facile scheme in which 
hereditary strains are computed using a recursive relation suggested by Henriksen. Composite 
response is studied using two competing micro-models, viz. a simplified Square Cell Model 
(SSCM) and a Finite Element based self-consistent Cylindrical Model (FECM). The algorithm is 
developed assuming that the material response computations are carried out in a module attached 
to a general purpose finite element program used for composite structural analysis. It is shown 
that the SSCM as used in investigations of material nonlinearity can involve significant errors in 
the prediction of transverse Young’s modulus and shear modulus. The errors in the elastic strains 
thus predicted are of the same order of magnitude as the creep strains accruing due to visco- 
elasticity. The FECM on the other hand does appear to perform better both in the prediction of 
elastic constants and the study of creep response. 



Nonlinear Visco-elastic Analysis of Composites 
via Micro-mechanical Models 


1. Introduction 


Highly flexible structural components are beginning to receive serious attention because of the 
increasing prospects of their use in aerospace vehicles. Such structures are fashioned out of materials 
which are liable to have time-dependent properties. In general, the response may be characterized as 
nonlinearly viscoelastic (NVE). 

There have been several attempts to a develop theory that would describe the NVE response of 
anisotropic or generally orthotropic composites. Significant contributions in this are due to by Schapery 
and his co-workers (Lou and Schapery, 1971; Schapery, 1974; Schapery and Sicking, 1995). The analysis 
requires several parameters that need to backed out from tests. To implement such an approach tests on 
lamina with a certain fiber volume fraction must be performed to identify the parameters that characterize 
its behavior. Such results are, however, applicable only to that material system with that fiber volume 
fraction. The variation of properties of the composite with change in fiber volume fraction (V f ) may be 
necessary for design optimization or simply to take into account the differences in V f between test 
specimens and prototypes, which may result due to variations in manufacturing process. 

On the other hand, there is definitely a higher level of confidence in our understanding of the 
NVE response of isotropic materials. A single integral law was first proposed by Schapery(1969) for 
uniaxial deformation in terms of instantaneous and creep compliance. The nonlinear response was 
inducted by letting the four parameters of the formulation ( go, gi, g2 and a CT ) be stress-dependent . This 
was subsequently generalized to account for 3-D states of stress, using the same parameters and letting 
them vary in the with the effective stress ( of the J 2 -plasticity theory) in the place of uniaxial stress in the 
same manner. Because of its simplicity and ability to describe the response of the material under complex 
loading history, it has been used by a number of investigators (e.g. Henriksen, 1984, Lai and Bakker, 
1996). 


In polymer-based composite, it is the matrix that exhibits significant nonlinearities and time- 
dependent behavior and polymers may be treated as isotropic. Given this, micro-mechanical models 
which treat fiber and matrix as distinct entities offer themselves as powerful tools for predicting the 
behavior of the composites, given the properties of the constituents. These offer a simpler and a more 
scientific approach to the problem in that only the basic constituents, viz. the fiber and the matrix need to 
be tested. TTie effect of varying fiber volume fractions and the influence of 3-D states of stress not 
captured in the tests are automatically accounted for by computer simulations. This obviates the need for 
postulating a priori multi-axial constitutive relationships of composites. Also the initiation of localized 
failure due to excessive strain can be accounted for in the micro-model with relative ease - something 
which orthotropic visco-elastic theories would be unable to do. The approach lends itself to a modular 
approach to material response, e.g. one can use different models such as plasticity, nonlinear visco- 
elasticity, visco-plasticity, etc. It is clear therefore that the micro-mechanics offers a scientific approach as 
against empirically based one to solve a variety of complex problems. However, the selected micro-model 
must be sufficiently simple and accurate for its widespread use in dealing with problems encountered in 
engineering practice. 
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Schaffer (1980) and Schaffer and Adams (1981) studied the micro-mechanical response using a 
square representative volume element (RVE) inside which a circular fiber is embedded. A detailed finite 
element analysis was performed and the analysis procedure required memory of all the previous stress 
history at any increment of loading. Consideration of normal and shear stresses running in the fiber 
direction would require 3-D elements in the place of 2-D plane strain elements used by Schaffer and 
Adams. Extension of such an analysis to a laminate under general loading conditions would involve 
installing such an analytical model at every integration point which would prove prohibitively expensive 
in terms of man/machine time and effort. 

The use of finite elements is not an insurance against poor accuracy of the results obtained. The 
postulated boundary conditions for the model do play a part in the solution obtained. Consider for 
example a model with square RVE; the loads across the fiber are usually applied by prescribing uniform 
displacements. In this case only an upper bound to the transverse modulus (E T ) would be obtained, 
however refined the mesh employed, the result converging to a value higher than that predicted by models 
where such restrictive boundary conditions are avoided. If on the other hand, the stress is prescribed, a 
lower bound is obtained. It appears best not to impose stresses or strains directly on the model, but to 
embed the model in a medium with properties which are the same or as close as possible to the 
homogenized properties of the composite itself and prescribe boundary conditions in the far-field. This 
leads us to the concept of self consistent cylindrical model first proposed by Hashin and Rosen (1964). 

A very simplified micro-model in the family of square cell models is that due to Aboudi (1980) 
referred to in literature as Aboudi Method of Cells (AMC). A simplified square cell model has been 
formulated by Sridharan and Jadhav (SSCM) which is equivalent to AMC and differs from it only in the 
manner of development and presentation (Sridharan and Jadhav, 1997) - the latter being employing only 
elementary mechanics arguments. However, because of the restrictive nature of simplifications involved, 
the stress distribution predicted by the model (AMC/SSCM) can be poor in accuracy (Sridharan and 
Jadhav, 1997) and this can lead to serious discrepancies in the prediction of nonlinear response of 
composite materials.(Jadhav and Sridharan, 1999). The conventional approach to minimizing the 
inaccuracy is to fit the parameters of the nonlinear response of the matrix such that the predicted 
composite response is in agreement with an experimentally determined key response, typically in shear ( 
Arenburg, 1988). Thus the approach tends to become somewhat empirical and loses its universal 
applicability. Starting with properties of “neat” resin and fibers, Jadhav and Sridharan (2002,2003) 
employed SSCM and a finite element based (self-consistent) cylindrical model (FECM) respectively for 
the analysis of carbon-epoxy laminae, laminates under off-axis loading and cylindrical shells under 
compression. Thermal residual stresses due to curing, matrix nonlinearity and cracking were all duly 
considered in developing the models. The results produced by SSCM were generally seen to be poor. The 
FECM, on the other hand fared much better in comparison to experimental results. The SSCM model, 
however still continues to be employed apparently because of its relative simplicity. Recently it has been 
used to study the nonlinear visco-elastic response of fiber reinforced plastics (FRP) by Muliana and Haj- 
Ali (2004). In the present work, we present the results of the SSCM along with those of FECM 


The self consistent cylindrical model as developed in the present study consists of an interior 
fiber element, which is surrounded by an annular (ring) matrix element which in turn is surrounded by an 
annular composite element of arbitrary width in the radial direction. The axial normal and shearing strains 
are taken as constant in the axial direction, which renders the problem two dimensional. The model is 
designated FECM: Finite Element based Cylindrical Model. As compared to the version presented earlier 
(Jadhav, 2000; Jadhav and Sridharan, 2002,2003), the model has been further simplified without a 
detectable loss of accuracy, making it more viable in the nonlinear finite element analysis of composite 
structures. 
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Fig. 1 RVE in the form of a circular fiber embedded in the matrix. 
Under transverse normal stresses 
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In the present work, both SSCM and FECM are extended to the study of nonlinear visco-elastic analysis. 
A significant feature of the micro-model analysis is that it employs the use of recursive relationships for 
the computation of hereditary effects (Henriksen, 1984) in the matrix and is developed as a material 
module for incorporation into a nonlinear finite element package. We note similar developments have 
been reported by Muliana and Haj-Ali in the context of SSCM ( Muliana and Haj-Ali, 2004 ). 

Essentia] Ingredients of Analysis : 

The essential ingredients of structural analysis scheme developed here are : 

(i) Constitutive relationship of the materials involved; the fiber may be treated as linearly elastic and the 
matrix as isotropic and nonlinearly viscoelastic. 

(ii) Micro-mechanical model used to transition from the material to the composite response. 

(iii) A geometrically nonlinear finite element scheme which incorporates the above mentioned features 
and is tailored to the study of long term structural behavior . 

The foregoing aspects will be first described in the sequel. Next the nexus between the material module 
and the global finite structural analysis would be described. 
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2. Nonlinear Visco-elasticity Model 


2.1 Behavior under Uniaxial Stress: 

The nonlinear viscoelastic behavior under uniaxial stress is simply described by an equation 
proposed by Schapery (1974) and used subsequently by a number of investigators (see for example 
Schaffer and Adams, 1981, Henriksen, 1984, Aboudi 1990, Sung Yi et. al , 1998). It is stated in the form: 

m = g'„ D a a(t) + gj (1) 

where: 

t is the time at which the total (kinematic) strain £ is given, 
a is the corresponding stress; 

D 0 is the initial value of the compliance at t = 0; 

go', g/, g 2 ‘ are the material parameters which dependent on the current stress a(t) and enhance the 
instantaneous and transient ( hereditary) strain contributions at time t; 

\\> is the reduced time given by: 

*<‘> = it <2) 

where a„ is a time scale shift factor, dependent on the stress a(s). Note that the superscripts t and s are 
introduced to indicate the current time t and the earlier times, s inside the integral. 


D c is the transient creep compliance which can be expressed in the form of a power law or in the form of a 
Prony Series: 


= D c if>" = S Z>, [1 - e ( - l ' r) ] 

r 


( 3 ) 


For a linearly viscoelastic material, the parameters go* = g/ = g 2 s = a c = 1 . The variation of the parameters 
go, gi , g 2 , and a<j with stress can be determined by creep and recovery tests. 


The treatment of nonlinear visco-elasticity follows Henriksen (1984) which is summarized in the sequel. 
In an incremental computation in reduced time domain, Eq.l can be expressed in an operational form : 

s' = F(a ‘ ) = Dj o' + E' ( 4 ) 

where the instantaneous compliance D\ is given by : 
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d) = g' 0 D 0 + g\g ' 2 Z£> r (i-r/ ) 


(5) 


where T* is a relaxation coefficient in the creep compliance series, defined as 

Ar “ A r A'F' 

where the reduced time increment is given by : 
and E‘ is the hereditary strain contribution given by : 
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where the q term is given by : 
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Assuming that gj <7 varies linearly over the current step, a recursive relationship (Henriksen, 1984) 
between the q’s at t and t-At can be written in the form: 


q r = 


t-At -XfA'V 1 

q r e r 


It _I „t-At 
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( 10 ) 


A detailed derivation is available (Sridharan, 2004). 

2.2 Isotropic material under multi-axial stress 

For isotropic materials, it may be assumed no shear-normal coupling develops as a result of 
multiaxial loading and the following nonlinear visco-elastic stress- strain relationship is proposed : 

{*'} = (H) 


where j^j 



is a vector of the algebraic difference between the kinematic strains 


w 


and thermal strains 
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( 12 ) 


{s‘J = {{e[-0') {e' 2 - 0 ') {4-0') e 4 e 5 e 6 } 


Here the single subscripted e’s stand for engineering strain components and are related to their tensorial 
counterparts as in: 


{ «1 e 2 e 3 

e 4 e 5 
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The expression for a certain E t 

takes the same form as E in eq.(8) where a is replaced by ctj (eq. (13.a) 


and q‘ r by 4 ', (eq. 10). Thus we have : 
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(15) 


where G\ = g' 2 erf 

[N] in eq. (1 1) is a 6 x 6 matrix defined in terms of Poisson’s ratio 
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In multi-axial stress, the 

parameters go , gi , 

g 2 , and a,, 

are taken as dependent upon the equivalent stress 


given by : 
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Inversion of eq.(l 1) gives a constitutive relationship 


< ,) 1 ( r ,1 f 

} =^7U M ‘Jn - 

where [M] = [N]' 1 . Alternatively, 

m m 

where [5; ] is the instantaneous material stiffness matri 
multiaxial case given respectively by 
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3. Micro-mechanical Models 


The purpose of the micro-models is to obtain the stress-strain state of the matrix and fiber given 
the stress-strain state of the composite and provide the stiffness of the composite, at the end of each load 
increment. In the present study, two competing micro-mechanical models are studied. They are briefly 
described in the sequel. 

3.1 Simple square cells based model, SSCM. 

It is assumed that fibers in the form of square bars are arranged in the form of doubly periodic 
array in the matrix. The representative volume element (RVE) for this model is a square cell, comprising 
of one fiber sub-cell and three matrix sub-cells (see Figure 1). X 2 and X 3 axes are lines of symmetry of the 
configuration and fibers runs along the X! axes. X 2 axis lies on the plane of the lamina and X 3 is the 
transverse axis. The orientation of the fiber axis to the reference axes of the lamina is supposed to be 
given. 


The fiber cell is a square of width, w f = ^v~ f (v f = fiber volume fraction); the two matrix cells 

contiguous to the fiber cell, have dimensions of w f x w m each while the third is a square of side w m , 

where w m = \- w f . The stresses and strains in the homogenized composite are equal to the volume 

averages of the respective quantities over these four sub-cells. The stresses and strains of the sub-cells are 
connected through a set of equilibrium and compatibility conditions (Sridharan and Jadhav, 1997; Jadhav 
and Sridharan, 2002, 2003) 

The equilibrium equations are : 
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In the foregoing equations, a bar represents an externally applied stress in the model or stress in the 
composite material at an integration point. The superscript (m,n) denote a certain cell in the model 
(Fig-2). 

The compatibility conditions are: 
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(2 1 .a-f) 

where a bar represents an externally applied strain in the model or a strain computed at an integration 
point in the composite material. 

To these equations, we must add 24 equations giving the 6x6 stress-strain relations in each cell. They 
completely describe the behavior of the model, viz. the stiffness and compliance of the model or the 
internal stresses under a set of external actions. Computations are carried out in terms of reduced 

unknowns, 10 in number : {x} T . 
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Fig. 2 Configuration of a SSCM 
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3.1.1 Determination of Elastic Constants : 

In order to determine the engineering elastic constants, unit values of 
, (T 22 ? ^33 5 23 >^"l3 ’^"l2j successively and determine the corresponding £y . This 

gives the compliance matrix from which the engineering constants may be found in the customary 
manner. 

It is important to note that the material has 6 independent constants, possessing as it does, tetragonal (or 
simply square) symmetry, in contrast to the anticipated transversely isotropic characterization involving 5 
independent constants. This is important and will be discussed in some detail in a later section. 

3.2 Finite Element based Cylindrical Model (FECM) : 

The RVE in this case is a composite cylinder (Fig.2) consisting of three concentric cylindrical elements. 
The outermost is the homogenous composite medium, the middle one is the matrix and the innermost is 
the fiber. Only the interior matrix and fiber elements constitute the model. The stresses and strains are 
volume averaged only over the two interior elements. The stresses however are applied at the boundary of 
the composite element. The composite stiffnesses are taken to be those that obtain at the previous 
increment or an estimate based on simpler models. 

The strain-displacement relations in polar coordinates for a three dimensional continuum are employed. 
The displacement functions are in the form of selected harmonics in the circumferential (0 -) direction 
and p-version type polynomials in the radial direction. The strains associated with z-direction are taken as 
constant. The fiber element is also treated to be annular (ring) element with a minute hole at the center. 

In the elastic range, the harmonic terms that are needed can be identified for any given loading 
case - either a sine or cosine harmonic for a certain displacement component. In the case of combined 
loading, the sine and cosine harmonics must be taken together. Since the displacements and stresses are 
not axisymmetric, the nonlinearity of the stress-strain relationship results in a non-uniform (non- 
axisymmetric) stiffness variation in the circumferential direction. This means that additional harmonics 
and harmonic coupling could come into play. However as seen in previous studies, (Jadhav and 
Sridharan, 1999,2002 ) the first order the effects of stiffness variation can be captured by accounting 
forthe variation of material stiffness in the circumferential ( 6 -) direction without augmenting the 
harmonic terms. A brief summary of the formulation is given in the sequel. 
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3.2. 1 Strain-displacement Relations: 
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where u, v and w are displacement components in the axial ( Z-direction , Xj), tangential and outward 
normal directions respectively; e’s are normal strains and y’s are shear strains in the appropriate senses 
indicated. 


3.2.2 Shape functions 

u = {zu o + U c u ^(r) cos(0) + < ^(r)sin(^)} 

v = k, ^, (r)cos(26>) + v'. (f> i (r ) sin(2^>] (23.a-c) 

w = [ w o, M r ) + w 2 i <P,(r)cos(20) + w 2 c , (j> i (r ) sin(2(9)] 

The shape functions for v and w enclosed within square brackets pertain to the action in the X 2 -X 3 plane 
and those for u enclosed by the braces are associated with the stress/strain associated with the Z-axis. 
However the u<, term in eq.(21.a) couples with the w 0 terms which describe the transverse action because 
of Poisson effect. 

It is seen that the shape functions are abridged from those presented previously (Sridharan and Jadhav 
,1997) and constitute the simplest version of the method. It will be seen that accuracy is not compromised 
as judged from the numerical results. 


3.2.3 Determination of Elastic Constants : 

In order to determine the engineering elastic constants, unit values of 

i , a 22 i (T33 , ^23 5 <3”i3 1 0"i2 } successively in the far-field. The volume averages of both 

stresses and strains are computed considering only the fiber and matrix elements. These give the 
necessary equations for obtaining the compliance matrix from which the engineering constants may be 
found in the customary manner. This model characterizes the material as transversely isotropic in terms of 
five independent constants. 
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3.2.4 Stress-strain relations as functions of 6. 


In an analysis with material nonlinearity, the material stiffness/compliance parameters are dependent on 
the current stress level and therefore vary spatially. In the present study, the instantaneous compliance 
parameter D ' of the matrix is dependent on g 0 ,g l , g 2 and T r which in turn depend on the effective 

stress at any location. It is convenient in computations to express Sj =1 / D\ in terms of a Fourier series in 

9. The variations of the parameter are not too localized and it is found a small number of terms should 
suffice. Keeping in view the couplings liable to occur in setting up of the equilibrium equations, we may 
express the variation of S] for any given r in the form: 


Sj = S 0 + ^ S c n cos( nO ) + ^ S* sin( nG) 

l l 

The strains are expressed in polar coordinates in the form: 

{*} = [*]{< 7 } 


(24) 


(25) 


where {q} are the degrees of freedom for unit stress applied in a certain direction. [B] matrix can be 
expressed in the form: 


[5] = [B 0 ] + [ff i ]cos( 0 ) + 1 ]sin( 9 ) 

+ [fi I ]cos( 20) + [»| ]sin( 20) 

The corresponding stresses are expressed as in eq.(17) in terms of instantaneous stiffnesses 

4 4 

{cr} = <Jq + cos(«#) + ^^ cr s n sin(«f?) 


(26) 


(27) 


Equilibrium equations take the form : 




r= RdO d z 


(28) 


where {p} are the boundary forces applied at r = R ( radius of the composite element) transformed to the 
polar system and {<(>} are the corresponding shape functions. 


The instantaneous stiffnesses of the model are computed in the same manner as the initial elastic 
constants by the application of unit stresses <Tj in turn at the boundary of the model, (setting the time 
variation to zero). Only difference is the variation of material stiffness point to point in the matrix. 
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Circumferential variation is handled using eq.(24) . The strains averaged over the interior two elements 
(matrix and fiber) yield corresponding instantaneous compliance terms. 

3.3 Determination of Hereditary Strains and in the Micro-models : 

In this section we outline the procedure for finding the hereditary strains and current stresses for the 
composite. (In fact the two depend on each other.) We shall assume that the hereditary strain distribution 
in the matrix are known at any given time. 

3.3.1 SSCM : 

Equations (14), (15) and (19.b) give the hereditary strain components {E ! } in terms of the current 
stresses at time t, and the stresses at time t - At. In SSCM, these strains in the three matrix cells will, in 
general assume, differing values. This will violate the compatibility conditions given by eq.l9(a-f). An 
auxiliary strain system {As a } must then be imposed so that 

(i) the compatibility conditions are satisfied and 

(ii) no additional stresses are imposed on the model. 

The current stresses in each cell take the form: 

V'} = [S}]|E'- A, +As-{P} + As a } (29) 

where {s , At } and {Ae} are the values of the strains in a certain cell as found at the end of the last 
increment and estimates of incremental strains at the current increment. [SJ] is the instantaneous 
stiffness. The unknowns {Ae a } are obtained invoking the governing equations of the model. The corrected 

hereditary strains are obtained as {E*} - {As a } for each cell. These strains are volume averaged to 

obtain the hereditary strains for the composite. The composite stresses are obtained from the current 
composite kinematic strain and the corrected composite hereditary strain as follows: 

k) = [ C S/ ] + Ae c - }} (30) 


where the subscript c indicates that the values pertain to the composite as a whole. 
3.3.2 FECM 


In this case, the hereditary strains in the matrix will vary along the circumferential as well as in the radial 
directions. The integration stations are selected to coincide with Gaussian integration points in the radial 
direction and at specified intervals (of say 15°) in the circumferential direction and the hereditary strain 
components are computed at each of the integration points using Eq.(14),(15) and (19.a). These strains are 
associated with additional displacements which must be so determined as to leave stresses in the 
composite unchanged. 


Let |CT 1 1 denote the stresses at time t at any point. Writing this in a form similar to Eq.(19) by including 
the effect of corrective displacements, w'e have, 
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\r'\ = + 


( 31 ) 


The basic equation of equilibrium is given by: 

tp) - JWW rdrdddz = 0 (32) 

vol 


where {p} is the vector of external forces, obtained from the current composite stresses. These external 
stresses which are in the global system are transformed to the polar system. The second term gives the 
internal nodal forces. Eq.(32) gives {Aq} from which the corrected hereditary strains 

{ {E 1 } - {B]{Aq }} are obtained at any given point. Note that in performing these calculations only the 

fiber and the matrix element alone need be considered and the averaged composite stresses at time t ( or 
the best estimate thereof) are deemed to be acting on the boundary of the matrix element. Volume 
average of the corrected hereditary strains over the two elements gives the corresponding composite value 

{E ‘ } . Equation (31) gives the stress distribution in the interior elements at time t. 


3.4 Nexus between Global Finite Element Program ( GFEP) and the material module 

The micro-models described in the foregoing are generally used as material modules appended to a 
general purpose commercial program ( called Global finite element program : GFEP, here) in order to 
solve large scale problems The interaction between GFEP and the material module is through an 
integration point at a certain location of the structure. For simplicity we assume that such integration 
points are located in each lamina represented by the micro-model across the thickness of a laminate. All 
calculations pertaining to visco-elasticity are done in the material module only and GFEP performs 
structural calculations on the homogenized material. 

For an iteration within an increment of load/time, the GFEP makes available to the micro-model the 
following quantities: The total strain components at the end of the last increment , i.e at t -At, and an 
estimate of the incremental strain components for the increment. Thus the current values of the total 

kinematic strains are Vc ^ +A£- c } .These must be kept unchanged in the material module and any 

values returned to the GFEP should be based on these strains. (In the last section, hereditary strain 
distribution as well as the internal stress distribution were corrected keeping the composite stress 
constant, but the purpose of this was just to determine the composite hereditary strains to be used in the 
module in the next incremen.) Thus the stresses returned to the GFEP must be computed ensuring the net 
changes in strain contributed by {Ae a } in SSCM and the corrective strains given by {Aq} using eq.(27) in 
FECM must be zero, maintained at values returned by GFEP, 

This is achieved in FECM by simply deactivating the degrees of freedom that are active on the boundary 
between the matrix element and the composite element. Equation (32) is set up by assembling the 
stiffness matrices for the interior two elements and the load vector involving the hereditary strains in the 
matrix. The degrees of freedom {Aq} active at the outer periphery of the matrix element are set to zero. 
Considering this with eq.(32), this means that the load vector {p} as well as the outer matrix element go 
out of action. The solution process consists of solving the equation : 
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0 


(33) 


J [B] T {cj}dv = 

vol 

for Aq with eq.(31) giving a, with boundary conditions mentioned already. The volume averaged 
stresses over the interior two elements are supplied to the GFEP. 

A step by siep procedure for computation is outlined in the sequel to complement the foregoing 
description of the models. 


3.5 Step by Step Procedure : 


The solution is sought incrementally, with the time at the beginning of an increment denoted by t - At 
and At being duration of the increment. 


The following quantities are stored as the values of state variables of the model : 

(i) and G\ ^ for each of the matrix locations which occur in the equation (14) for the 

. ,. r’t — At 

hereditary strains bj , 

(ii) the instantaneous parameter Dj~ M for the matrix locations as well as the composite stiffness 

[*«/"*] 

(iii) the effective hereditary strains for the composite ^ j 

The “ matrix locations" referred to above are the three matrix cells for SSCM and for FECM, these are 
mo x ns stations where mo is the chosen number of Gaussian points for integration in the radial direction 
and ns is the chosen number of equi-spaced integration stations in the 9-direction. 

The following steps are involved in the analysis: 


1. Using the matrix and composite material properties available at t- At, and the fiber elastic 
properties determine the composite compliances by applying unit values of stress component one 
at a time and determining the volume averaged strains. Transform the compliances to 
stiffhesses.(For SSCM, it is found convenient to apply unit strain in the z-direction instead of 
averaged stress and the procedure is modified accordingly.) 

2. Transform the global strains given by GFEP to the local system (the axes of material symmetry). 

3. Update the composite stresses from {d :t } = [ C »S/ + As c — E* c ^ j 

(This involves an approximation which can be corrected by iterating.) 

4. Determine the stresses in the fiber and matrix locations in the model using the micro-model 
equations. 

5. For each matrix location: 
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(i) compute g‘ 0 ,g[ ,g' 2 ,a' a ,A<//' ,r ; , ( i = 1,..N) , the instantaneous compliance parameter D \ ; 

(ii) update q’s using recursive relationships ( Eq. 10)andG’s finding q' rlj and Gy and 

compute the hereditary strains E f at all matrix locations, eq.(14). 

6. From the instantaneous stiffness of the matrix and fiber, the corresponding instantaneous 

stiffness of the composite [ C S/] is determined. This is transformed to the global coordinate 
system and returned to the GFEP. 

7. Computation of averaged hereditary strain \ e 1 c j : The stresses are taken in the form given by Eq. 

(29) or (31) and the set of corrective strains or corrective degrees of freedom are so computed as 
to satisfy the micro-model equations (32)and contribute no additional averaged stresses. This 
gives the corrected hereditary strain distribution (vide section 3.3) which averaged over the 
volume gives the composite hereditary strain .(This is for use in the next increment). 

8. Find the stresses at time t, cf/ to be returned to GFEP. As per Art. 3.4 , a set of corrective 

kinematic strains are obtained which satisfy compatibility and average to zero over the micro- 
model. These are employed in Eq(29) for SSCM and Eq.(31) for FECM, to give the updated 
stresses. These are volume averaged , transformed to the global system and returned to GFEP. 


These calculations are performed every iteration till convergence is achieved for the increment. 


3.6 Determination of Engineering Elastic constants by SSCM - revisited 


Before the more involved nonlinear visco-elastic analysis is performed, it is incumbent on us to verify the 
accuracy of the models in the prediction of linear elastic behavior - something that is too often taken for 
granted. This will help isolate errors due to the micro-model employed from those due to inherent 
limitations, if any, of the Schapery theory. 

Versions of SSCM : 

As far as determination of elastic constants of the composite material is concerned, it is helpful to 
distinguish , in the interests of clarity, between two versions of SSCM, one used in finding stresses 
carried by the matrix and nonlinear analysis in general (Version I) and the other used exclusively for the 
determination of elastic constants (Version II). 

Version I : As already mentioned, the use of the governing equations (20.a-n) and (21.a-f) , produces a 
material with six independent constants and the transverse shear modulus (G T ) is not related to Transverse 
Young’s modulus (E T ) and the transverse Poisson ratio (v T ). The RVE which is square in its outline is 
assumed to be so oriented as to have its edges parallel to the principal axes of the material. 


Version II : The last assumption is relaxed in this version, i.e., the RVE’s representing a small 

neighborhood of the material are all randomly oriented with respect to the principal axes of the material. 
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A typical RVE would make an angle 0 with respect to the X 2 -X 3 axes (Fig. 5 ) where 0 may vary from 0 
to it. Thus the elastic constants in the X 2 -X 3 plane must be averaged by appropriate integration from 0 to 
p. The results of this operation are give the [C’] the corrected material stiffness matrix which 
characterizes the material as transversely isotropic : 


C'li 

= 

C \\ 


C'n 

= 

C'n = Ci 2 


C 22 

= 

C'33 = ^ C22 

H — C 73 H — C44 
4 23 2 44 

£”23 


T C 22 + t £'23 “ 
4 4 

1 C 
2 44 

C -44 

= 

2 ( C 22 -C 23 ) 

= 2^ 22 ~ ^23 ) 

£”55 

_ 

C 66 = C 55 



The engineering elastic constants can be obtained from well known relations of elasticity. While it makes 
sense thus to average the stiffnesses, the stresses cannot be so averaged meaningfully and as such this 
version is limited to finding the elastic constants. It is important to bear in mind this dichotomy, as the 
accuracy of Version I cannot be established on the basis of accuracy of Version II. 
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4. Numerical Examples 


4.1 SSCM vs FECM : A Comparison of Elastic Constants 

Two typical composite materials are considered, viz.. Glass-epoxy and Carbon Epoxy. The glass is 
isotropic whereas carbon is transversely isotropic. The properties are listed in Tables 1 and 2. 


Fig.6(a-b) show the variation of transverse Young’s Modulus (E T ) and shear modulus (G T ) respectively 
of Glass-Epoxy for fiber volume fraction ranging from 0.4 to 0.8. Predictions of both versions of SSCM 
as well those of FECM are shown in the figures. Fig.6c. shows the variation of axial shear modulus (G A ) 
with V f . Note the versions I and II of SSCM produce identical results as far as axial properties are 
considered. The results of E T and G A are compared with the elasticity solution by Pickett( Aboudi, 1991) 
in Fig.6.a and 6.c respectively. 

Similar results are displayed in Fig 7(a-c) for carbon-epoxy. These results are compared with the 
experimental results given by Kriz and Stinchcomb ( Aboudi, 1991). 


Consider first the predictions for E T of glass epoxy. (Fig.6.a). It is seen that there is a significant 
discrepancy in the prediction of E T between the two versions of SSCM. The version II is in very good 
agreement with FECM as also with the elasticity solution. The discrepancy in the Version I prediction is 
of the order of 15-20% with Version II as the datum. Such a serious difference in the elastic problem 
raises a question of the quality of results produced by SSCM in a nonlinear analysis which is based on 
Version I. More significant differences are noticed in the prediction of G T which are of the order of 30%. 


The SSCM and FECM are in good agreement in the prediction of axial shear modulus, G A and these in 
turn are in agreement with the elasticity solution (Fig. 6.c). The use of a single variable u, in contrast to 
earlier formulation (Sridharan and Jadhav, 1999) involving all the three displacement components in the 
analysis of axial shear response has produced results of requisite accuracy. 

Next, consider the case of carbon epoxy. The discrepancies between the two versions for carbon-epoxy 
are considerably less significant and are about 5% (Fig.7.a-b).This is because the difference in the 
modulus of the matrix and the transverse modulus of carbon epoxy is much smaller than in the case of 
glass epoxy. As before the FECM and Version II are close to each other. The experimental results for 
transverse shear modulus clearly favor the latter two, but that for Et (Fig. 7.a) tends to favor Version I. 

Taken all in all, FECM produces consistently accurate results for both Glass-Epoxy and Carbon Epoxy. 
The results of Version I must be viewed with caution for Glass-epoxy. 
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Table 1 : Elastic constants of Glass fiber and Epoxy matrix, case (i) 


E (GPa) v 

Glass 

68.94 

0.2 

Epoxy 

3.42 

0.34 


Table 2 : Elastic constants of for Modmor Carbon fibers and epoxy matrix. Case (ii) 


Phase 

E a 

(GPa) 

Va 

E t 

(GPa) 


G a 

(GPa) 

Carbon 

(Transversely 

Isotropic) 

232 

0.279 

15.0 

0.49 

24.0 

Epoxy 
( Isotropic) 

5.35 

0.354 ! 

5.35 

0.354 

1.975 
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-Epoxy with fiber volume fraction 





Transverse Shear Modulus, Gp, GPa 





0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 

Fiber Volume Fraction, V f 


Fig. 6(c) Variation of Axial Shear modulus, G T of Glass -Epoxy with fiber volume fraction 
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Transven 



Fiber Volume Fraction, V f 


Fig. 7(a) Variation of Transverse modulus, E T of Carbon-Epoxy with fiber volume fraction 




Transverse Shear Modulus, Gp, GPa 



Fiber Volume Fraction, V f 


Fig. 7(b) Variation of Transverse Shear modulus, Gj of Carbon-Epoxy 
with fiber volume fraction 
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Fiber Volume Fi 


Fig. 7(c) Variation of Axial Shear modulus, G T of 
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4.2 Creep of Glass-epoxy Composite 

As mentioned earlier Schaffer (1980) studied both experimentally and analytically the time-dependent 
behavior of composites. He considered, among others, Hercules 3501-6 epoxy resin which exhibits visco- 
elastic effects reinforced by S2 glass fibers. The following material properties are available: 

Glass is considered isotropic with E= 86.9 GPa, and v = 0.22. 

The initial value of the E of the epoxy = 4.14 GPa and v = 0.34. 

The fiber volume ratio, V f = 0.63. 

The stress-dependent nonlinear visco-elastic parameters have been measured and documented by 
Schaffer(1980) and as also by Aboudi (1991). The curve fitted expressions for g' 0 ,g' t ,g' 2 ,a' a are given 

in the appendix as also the creep compliance parameters ( D c and v, eq.(3) ) and the exponents ( Vs) and 
the coefficients (D r ’s) of the Prony series( Eq. 3). 

4.2.1 Case Study; 

Case (if : Creep of lamina under biaxial compression in the two transverse directions, given by 
a 2 =103.4 MPa and a 3 = 34.5 MPa. 

Case (if) : Creep of lamina under uniaxial compression at a stress level, a 2 , of 158.5 MPa. 

Results and Discussion - Case (i) ; 

Fig. 8 shows the following results which pertain to case(i). 

(1) Creep predictions given by Schaffer (1980) using a detailed finite element model which assumes the 
composite to consist of circular cylindrical fibers arranged in a doubly periodic array. The RVE cross 
section consists of a circle embedded in square region with the rest of the space filled by matrix. 

(2) Predictions based on the method of cells as reported by Aboudi(1991) 

(3) Predictions based on SSCM as detailed in this report.. 

(4) Predictions based on FECM 

It is seen that among the set of four results, Schaffer’s finite element analysis predicts the smallest values 
of strains throughout; in contrast the FECM predicts the highest values. The Schaffer’s analysis involves 
stipulating uniform displacement on the edges of the square, the boundaries of the representative area. 
Once a doubly symmetric fiber arrangement is assumed, specification of uniform displacements is 
inescapable, by virtue of the nature of symmetry involved. But as mentioned earlier, the specification of 
displacements in the analysis leads to an upper bound predictions of stiffness. (This is true of AMC and 
SSCM as well). In FECM no boundary conditions are prescribed directly on the model which is 
embedded in the surrounding composite medium. Thus Schaffer’s analysis results in smaller values of 
strain than that given by FECM at the instant of load application. However, the rate of increase of strains 
are not noticeably different in the predictions of the two models. 
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The results quoted by Aboudi(1991) and those of the present SSCM are quite close to each other. It is not 
clear why the two results do not coincide as they in essence employ the same set of equilibrium and 
compatibility conditions for the micro-model. The numerical treatment of the equations is, however, 
different. The present work is based on the use of recursive relationships in the computation of hereditary 
strains whereas in Aboudi’s analysis the entire past history of stress explicitly enters the calculations of 
hereditary strain at any moment in time. It is believed that the discrepancy is numerical in origin, e.g. a 
coarser time increment in Aboudi’s analysis. 

The differences between the square cell models and Schaffer’s model are the result of differences in the 
geometry. In the former, there are always two matrix cells in series in the resistance of transverse stresses 
which tend to make these models somewhat more compliant than Schaffer’s model. 

Case (ii) 

Fig. 9 plots the variation of (percent) transverse strain with time as obtained from tests and analyses, for 
case (ii). 

The following results are presented : 

1. and 2. : Experimental results reported by Schaffer(1980) 

3. FECM Results 

4. Results based Schafer’s model ( Schafer’s Prediction I) 

5. Schafer’s model with curing stresses accounted for ( Schafer’s Prediction II), and 

6. SSCM results 

It is seen that the two experimental results of nominally the same specimen and same test conditions are 
markedly different. The abrupt jump in the experimental results (II) at around 7 hours is somewhat 
suspicious. Perhaps the result designated as I is more representative of the actual behavior. 

Among the analytical results, Schafer’s model (Prediction I) gives the smallest value of strain at 
any time whereas the present FECM model gives the highest values. This trend is the same as for case (i). 
Accounting for thermal stresses due to curing increases the predicted values of strain by about 10%. It 
appears that the boundary conditions exaggerate this effect. The results of SSCM are once again higher 
than that given by Schafer’s model ( Prediction I). 

These results do not give any definitive indication regarding the relative accuracies of the models 
employed here. In view of the inaccuracies involved in Version I of the SSCM, the author’s previous 
experience with the micro-models (Jadhav and Sridharan, 2002,2003) and the intrinsic characteristics of 
the model themselves, it appears FECM is the more accurate model. 
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0 30 60 90 120 150 

t ( hours ) 

Fig.8. Creep of S2- glass -Hercules 3501-6 epoxy 
under Biaxial Compression, Case (i). 
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5. Concluding Remarks 


A general scheme for the analysis of nonlinearly visco-elastic composite materials which exploits the 
recursive relationships for the calculation of hereditary strains has been presented, together some 
preliminary results. The deficiencies of the simple square cell model are brought out by comparison with 
the self-consistent cylindrical model and the experimental results. There appears to be a broad qualitative 
agreement between the results of simpler analytical models developed here and the experimental results 
available in literature. More experimental results of greater consistency and less scatter are needed to 
validate the algorithm developed here for the prediction of nonlinear visco-elastic response of composites. 
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t ( hours) 

Fig. 9. A comparison of experimental results and the results given by analytical models 
for S2 Glass -Hercules 3501-6 epoxy under transverse compression. Case (ii). 






Appendix : Nonlinear Visco-elastic Material Parameters for Glass-Epoxy 

Given below are the curve-fitted expressions for the material parameters for the nonlinear visco-elastic 
material model , viz. go, gi, g 2 and a 0 . These are given in terms of the equivalent stress cr e . 

= 1.0 K<17); g 0 = 1.0 + 0.09 (<r e -17)/121 (<7, >17) 

g, = 1.0 (cr e < 17 ); g, = 1.0 + 2.38 (<r e - 17)/121 (<r e >17) 

g 2 = 1.0 (cr e <17 ); 

g 2 = 1.0-1.803(cr e -17)xl0" 2 +1.744 (<r e -17) 2 xl0" 4 -0.3847 (o- e -17>cl0- 6 (cr e >17) 

= 1.0 (cr e < 17 ); 

= 1.0-4.222 (< 7 e -17>cl0- 2 + 9.581 (a e -17) 2 xl0 -4 -0.1075(£7 e -17>d0' 4 (cr e >17 ) 

The creep compliance function( Eq. 3) takes the form : 

D c = 6.293x1 O' 6 ; n = 0.26 units involved are MPa-hours. 

The coefficients of Prony series D r and the corresponding X’s are given in Table. A1 


Table Al. Coefficients in the Prony Series 


r 

K 

D r 

i 

10° 

0.766139xl0' 5 

2 

10' 1 

0.405547x1 0' 5 

3 

10' 2 

0.1 11955x10" 

4 

10' 3 

0.16913x10" 

5 

10" 

0.39421xl0' 5 

6 

10' 5 

0.281654xl0' 2 

7 

10' 6 

-0.284248x10'' 

8 

10' 7 

0.407103x10' 

9 

10' 8 

0.701062xl0' 2 
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